Fit density (also referred to as CPUE) model with environmental predictors and use that to calculate weighted mean dissolved oxygen, temperature and depth of Baltic cod
# Load libraries, install if needed
library(tidyverse)
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(qwraps2)
library(wesanderson)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)
library(brms)
library(sdmTMBextra)
# To load entire cache in interactive r session, do:
# qwraps2::lazyload_cache_dir(path = "R/analysis/density_model_cache/html")
# Specify map ranges
ymin = 52; ymax = 60; xmin = 10; xmax = 24
map_data <- rnaturalearth::ne_countries(
scale = "medium",
returnclass = "sf", continent = "europe")
# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
sf::sf_use_s2(FALSE)
# https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data
swe_coast <- suppressWarnings(suppressMessages(
st_crop(map_data,
c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))
# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)
# Define plotting theme for main plot
theme_plot <- function(base_size = 11, base_family = "") {
theme_light(base_size = base_size, base_family = "") +
theme(
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.margin = margin(0, 0, 0, 0),
legend.box.margin = margin(-5, -5, -5, -5),
strip.text = element_text(size = 9, colour = 'gray10', margin = margin(b = 1, t = 1)),
strip.background = element_rect(fill = "grey95"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
}
# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 11, base_family = "") {
theme_light(base_size = base_size, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 7),
strip.text = element_text(size = 8, colour = 'gray10', margin = margin(b = 1, t = 1)),
strip.background = element_rect(fill = "gray95"),
legend.direction = "horizontal",
legend.margin = margin(1, 1, 1, 1),
legend.box.margin = margin(0, 0, 0, 0),
legend.key.height = unit(0.4, "line"),
legend.key.width = unit(2, "line"),
legend.spacing.x = unit(0.1, 'cm'),
legend.position = "bottom",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
}
# Make default base map plot
xmin2 <- 303000; xmax2 <- 916000; xrange <- xmax2 - xmin2
ymin2 <- 5980000; ymax2 <- 6450000; yrange <- ymax2 - ymin2
plot_map_raster <-
ggplot(swe_coast_proj) +
xlim(xmin2, xmax2) +
ylim(ymin2, ymax2) +
labs(x = "Longitude", y = "Latitude") +
geom_sf(size = 0.3) +
theme_plot() +
guides(colour = guide_colorbar(title.position = "top", title.hjust = 0.5),
fill = guide_colorbar(title.position = "top", title.hjust = 0.5))
d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/catch_q_1_4.csv")
# Calculate standardized variables
d <- d %>%
filter(year < 2020 & year > 1992 & quarter == 4) %>%
mutate(oxy_sc = oxy,
temp_sc = temp,
depth_sc = depth,
) %>%
mutate_at(c("oxy_sc", "temp_sc", "depth_sc"),
~(scale(.) %>% as.vector)) %>%
mutate(year = as.integer(year)) %>%
drop_na(oxy, depth, temp)
plot_map_raster +
geom_point(data = d, aes(X*1000, Y*1000))
pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
pred_grid <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")
pred_grid <- bind_rows(pred_grid1, pred_grid)
# Standardize data with respect to data grid:
pred_grid <- pred_grid %>%
mutate(year = as.integer(year)) %>%
filter(year %in% c(unique(d$year))) %>%
mutate(depth_sc = (depth - mean(d$depth))/sd(d$depth),
temp_sc = (temp - mean(d$temp))/sd(d$temp),
oxy_sc = (oxy - mean(d$oxy))/sd(d$oxy)) %>% # Need to scale these to the mean and sd in the data!
drop_na(oxy, depth, temp)
spde <- make_mesh(d, xy_cols = c("X", "Y"),
n_knots = 200,
type = "kmeans", seed = 42)
# Plot and save spde
png(file = "figures/supp/density/spde.png", units = "in", width = 6.5, height = 6.5, res = 300)
plot(spde)
dev.off()
# Depth spline + oxy spline
# Takes about 30 minutes
m <- sdmTMB(density ~ 0 + as.factor(year) + s(depth_sc) + s(oxy_sc) + s(temp_sc),
data = d,
mesh = spde,
family = tweedie(link = "log"),
spatiotemporal = "AR1",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_steps = 1))
#> Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
#> TMB was built with Matrix version 1.5.4
#> Current Matrix version is 1.5.3
#> Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
#> Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
#> glmmTMB was built with TMB version 1.9.3
#> Current TMB version is 1.9.4
#> Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
tidy(m, conf.int = TRUE)
sanity(m)
samps <- sdmTMBextra::predict_mle_mcmc(m, mcmc_iter = 201, mcmc_warmup = 200)
#> Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
#> https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
#> Warning: Examine the pairs() plot to diagnose sampling problems
mcmc_res <- residuals(m, type = "mle-mcmc", mcmc_samples = samps)
png(file = "figures/supp/density/qq.png", units = "in", width = 6.5, height = 6.5, res = 300)
qqnorm(mcmc_res);qqline(mcmc_res)
dev.off()
d$residuals <- mcmc_res
# Residuals on map
plot_map_raster +
geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = residuals), size = 0.5) +
scale_colour_gradient2() +
facet_wrap(~ year, ncol = 5) +
labs(color = "residuals") +
theme_facet_map() +
geom_sf(size = 0.1)
#> Warning: Removed 8 rows containing missing values (`geom_point()`).
ggsave("figures/supp/density/spatial_resid.png", width = 6.5, height = 8.5, dpi = 600)
#> Warning: Removed 8 rows containing missing values (`geom_point()`).
rho is ar_phi on
the -1 to 1 scale):tidy(m, effects = "ran_pars", conf.int = TRUE)
#> Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
#> TMB was built with Matrix version 1.5.4
#> Current Matrix version is 1.5.3
#> Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
#> Standard errors intentionally omitted because they have been calculated in log
#> space.
predict_mcod <- predict(m, newdata = pred_grid)
# Save prediction as csv (smaller than the model object...)
write.csv(predict_mcod, "output/predict_mcod_density.csv")
# Plot predicted density and random effects
plot_map_raster +
geom_raster(data = filter(predict_mcod, depth < 120 & depth > 10), aes(x = X * 1000, y = Y * 1000, fill = exp(est))) +
scale_fill_viridis_c(trans = "sqrt") +
facet_wrap(~ year, ncol = 5) +
labs(fill = expression(kg/km^2)) +
theme_facet_map() +
geom_sf(size = 0.1)
#> filter: removed 42,282 rows (16%), 216,297 rows remaining
ggsave("figures/supp/density/est_map.png", width = 6.5, height = 8.5, dpi = 600)
# Plot spatiotemporal random effect
plot_map_raster +
geom_raster(data = predict_mcod, aes(x = X * 1000, y = Y * 1000, fill = epsilon_st)) +
scale_fill_gradient2() +
facet_wrap(~ year, ncol = 5) +
theme_facet_map() +
geom_sf(size = 0.1)
ggsave("figures/supp/density/epsilon_st_map.png", width = 6.5, height = 8.5, dpi = 600)
# Plot spatial random effect
plot_map_raster +
geom_raster(data = filter(predict_mcod, year == 1999), aes(x = X * 1000, y = Y * 1000, fill = omega_s)) +
scale_fill_gradient2() +
geom_sf(size = 0.1)
#> filter: removed 249,002 rows (96%), 9,577 rows remaining
ggsave("figures/supp/density/omega_s_map.png", width = 6.5, height = 6.5, dpi = 600)
get_index_simsdiv_index_list <- list()
for(i in unique(pred_grid$sub_div)){
pred_grid_sub <- pred_grid %>% filter(sub_div == i & depth > 10 & depth < 110)
pred_sim_sub <- predict(m, newdata = pred_grid_sub, nsim = 500)
index_sim <- get_index_sims(pred_sim_sub,
area = rep(4*4, nrow(pred_sim_sub))) %>%
mutate(sub_div = i)
div_index_list[[i]] <- index_sim
}
#> filter: removed 221,103 rows (86%), 37,476 rows remaining
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> mutate: new variable 'sub_div' (double) with one unique value and 0% NA
#>
#> filter: removed 188,109 rows (73%), 70,470 rows remaining
#>
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> mutate: new variable 'sub_div' (double) with one unique value and 0% NA
#>
#> filter: removed 204,687 rows (79%), 53,892 rows remaining
#>
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> mutate: new variable 'sub_div' (double) with one unique value and 0% NA
#>
#> filter: removed 240,516 rows (93%), 18,063 rows remaining
#>
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> mutate: new variable 'sub_div' (double) with one unique value and 0% NA
#>
#> filter: removed 230,877 rows (89%), 27,702 rows remaining
#>
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
#> mutate: new variable 'sub_div' (double) with one unique value and 0% NA
div_index <- bind_rows(div_index_list)
# In addition to the subdivision indicies, we'll also group a few and also do the total
# Now get the total index by specifying the area of a grid cell
# Total prediction
preds_cod_sim <- predict(m, newdata = pred_grid, nsim = 500)
# Now calculate the biomass index
index_sim <- get_index_sims(preds_cod_sim, area = rep(4*4, nrow(preds_cod_sim)))
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
# Join indices
index_sim <- index_sim %>% mutate(sub_div = "Total")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
div_index_sim <- bind_rows(div_index %>% mutate(sub_div = as.character(sub_div)),
index_sim %>% mutate(sub_div = "Total"))
#> mutate: converted 'sub_div' from double to character (0 new NA)
#> mutate: no changes
# Now do the same but excluding the areas not sampled in the pred grid
ncells_sub <- filter(pred_grid, year == max(pred_grid$year) & depth > 10 & depth < 120) %>% nrow()
#> filter: removed 250,568 rows (97%), 8,011 rows remaining
pred_grid_sub <- filter(pred_grid, depth > 10 & depth < 120)
#> filter: removed 42,282 rows (16%), 216,297 rows remaining
pred_avg_sim_sub <- predict(m, newdata = pred_grid_sub, nsim = 500)
index_avg_sim_sub <- get_index_sims(pred_avg_sim_sub, area = rep(4*4, nrow(pred_avg_sim_sub)))
#> We generally recommend using `get_index(..., bias_correct = TRUE)`
#> rather than `get_index_sims()`.
# Combine and plot & compare
index_avg_sim_comp <- bind_rows(mutate(index_sim, area = "full"),
mutate(index_avg_sim_sub, area = "subset")) %>%
mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> mutate: new variable 'est_t' (double) with 54 unique values and 0% NA
#> new variable 'lwr_t' (double) with 54 unique values and 0% NA
#> new variable 'upr_t' (double) with 54 unique values and 0% NA
ggplot(index_avg_sim_comp, aes(year, est_t, ymin = lwr_t, ymax = upr_t, color = area)) +
geom_line() +
geom_ribbon(alpha = 0.4)
lab_size <- 1.7
div_index_sim <- div_index_sim %>% mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000)
ind_plot_sd <-
ggplot() +
scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
geom_line(data = filter(div_index_sim, !sub_div == "Total"),
aes(year, est_t, color = sub_div)) +
geom_ribbon(data = filter(div_index_sim, !sub_div == "Total"),
aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
guides(color = "none") +
labs(x = "Year", y = "Tonnes") +
theme_plot() +
guides(colour = "none",
fill = guide_legend(title.position = "top", title.hjust = 0.5)) +
theme(legend.position = c(0.6, 0.9),
legend.direction = "horizontal",
legend.key.height = unit(0.03, "cm"),
legend.key.width = unit(0.4, "cm"),
legend.text = element_text(size = 6),
legend.title = element_text(size = 7),
legend.background = element_blank(),
plot.margin = unit(c(0, 0, 0, 0), "cm"),
panel.background = element_blank())
ind_plot_tot <-
ggplot() +
geom_line(data = filter(div_index_sim, sub_div == "Total"),
aes(year, est_t)) +
geom_ribbon(data = filter(div_index_sim, sub_div == "Total"),
aes(year, ymin = lwr_t, ymax = upr_t, fill = "Total"), alpha = 0.4) +
scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
scale_fill_manual(values = "grey40") +
labs(x = "Year", y = "Tonnes", fill = "") +
theme_plot() +
theme(legend.position = c(0.88, 0.92),
legend.key.height = unit(0.03, "cm"),
legend.key.width = unit(0.4, "cm"),
legend.text = element_text(size = 6),
legend.title = element_text(size = 7),
legend.background = element_blank(),
plot.margin = unit(c(0, 0, 0, 0), "cm"))
cpue_map_94_18 <- plot_map_raster +
geom_raster(data = filter(predict_mcod, year %in% c("1994", "2018") & depth < 120 & depth > 10), aes(x = X * 1000, y = Y * 1000, fill = exp(est))) +
scale_fill_viridis_c(trans = "sqrt") +
facet_wrap(~ year, ncol = 2) +
labs(fill = expression(kg/km^2)) +
theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.5, "cm"),
legend.position = c(0.9, .085),
legend.direction = "horizontal",
legend.background = element_rect(fill = NA),
legend.text = element_text(size = 5),
legend.title = element_text(size = 7)) +
guides(fill = guide_colorbar(title.position = "top", title.hjust = 0.9)) +
geom_sf(size = 0.1) +
annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)
(ind_plot_sd | ind_plot_tot) / cpue_map_94_18 + plot_annotation(tag_levels = 'A')
ggsave("figures/density.pdf", width = 17, height = 15, units = "cm")
# All years cpue
plot_map_raster +
geom_raster(data = filter(predict_mcod, depth < 120 & depth > 10), aes(x = X * 1000, y = Y * 1000, fill = exp(est))) +
scale_fill_viridis_c(trans = "sqrt") +
facet_wrap(~ year, ncol = 5) +
labs(fill = expression(kg/km^2)) +
theme_facet_map() +
geom_sf(size = 0.1)
ggsave("figures/supp/density/all_years_density_pred.png", width = 6.5, height = 8.5, dpi = 600)
# Plot bathymetry
plot_map_raster +
geom_raster(data = filter(predict_mcod, year == "2006"), aes(x = X * 1000, y = Y * 1000, fill = depth)) +
scale_fill_viridis(option = "A", direction = -1) +
labs(fill = "Depth") +
theme(plot.margin = unit(c(0, 0, 0.5, 0), "cm"),
legend.key.height = unit(0.1, "cm"),
legend.key.width = unit(0.7, "cm"),
legend.position = c(0.8, .07),
legend.direction = "horizontal",
legend.background = element_rect(fill = NA)) +
geom_sf(size = 0.1) +
annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)
#> filter: removed 249,002 rows (96%), 9,577 rows remaining
lab_size <- 1.6
## Depth
# Plot bathymetry
bathy_plot <- plot_map_raster +
geom_raster(data = filter(predict_mcod, year == "2006"), aes(x = X * 1000, y = Y * 1000, fill = depth)) +
scale_fill_viridis(option = "A", direction = -1) +
labs(fill = "Depth", x = "", y = "") +
theme_plot(base_size = 11) +
theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.key.height = unit(0.1, "cm"),
legend.key.width = unit(0.35, "cm"),
legend.position = c(0.8, .085),
legend.direction = "horizontal",
legend.background = element_rect(fill = NA),
legend.text = element_text(size = 5),
legend.title = element_text(size = 6)) +
guides(fill = guide_colorbar(title.position = "top", title.hjust = 1)) +
geom_sf(size = 0.1) +
annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)
#> filter: removed 249,002 rows (96%), 9,577 rows remaining
# Now calculate quantiles of the biomass-weighted values
pal <- c(rev(brewer.pal(n = 8, name = "Dark2")), "gray20")
wm_depth <- predict_mcod %>%
group_by(year) %>%
summarise("Median" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.5)),
"1st" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.25)),
"3rd" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.75))) %>%
pivot_longer(cols = c("Median", "1st", "3rd"),
names_to = "series", values_to = "depth") %>%
group_by(series) %>% # standardize within for easy plotting
mutate(depth_ct = depth - mean(depth)) %>%
ungroup()
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> pivot_longer: reorganized (Median, 1st, 3rd) into (series, depth) [was 27x4, now 81x3]
#> group_by: one grouping variable (series)
#> mutate (grouped): new variable 'depth_ct' (double) with 79 unique values and 0% NA
#> ungroup: no grouping variables
plot_w_dep <- ggplot(wm_depth, aes(year, depth, color = series, group = series, fill = series)) +
stat_smooth(method = "gam", formula = y ~ s(x, k = 4), se = FALSE, size = 0.8) +
geom_point(size = 1.3, alpha = 0.8) +
scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
scale_color_manual(values = pal, name = "Quantiles") +
scale_fill_manual(values = pal) +
guides(fill = "none", color = "none") +
labs(y = "Depth [m]", x = "Year", color = "") +
scale_y_reverse() +
theme_plot(base_size = 11) +
theme(plot.margin = unit(c(0.9, 0, 0, 0), "cm"))
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
## Oxygen
# Plot oxygen in 2006
oxy_plot <- plot_map_raster +
geom_raster(data = filter(predict_mcod, year == "2006"), aes(x = X * 1000, y = Y * 1000, fill = oxy)) +
scale_fill_viridis() +
labs(fill = expression(paste("O" [2], " [ml/L]", sep = "")), x = "") +
theme_plot(base_size = 11) +
theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.key.height = unit(0.1, "cm"),
legend.key.width = unit(0.35, "cm"),
legend.position = c(0.8, .085),
legend.direction = "horizontal",
legend.background = element_rect(fill = NA),
legend.text = element_text(size = 5),
legend.title = element_text(size = 6)) +
guides(fill = guide_colorbar(title.position = "top", title.hjust = 1)) +
geom_sf(size = 0.1) +
annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)
#> filter: removed 249,002 rows (96%), 9,577 rows remaining
# Plot biomass-weighted oxygen
# Weighted total
wm_oxy <- predict_mcod %>%
group_by(year) %>%
summarise("Median" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.5)),
"1st" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.25)),
"3rd" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.75))) %>%
pivot_longer(cols = c("Median", "1st", "3rd"),
names_to = "series", values_to = "oxy") %>%
ungroup() %>%
group_by(series) %>% # standardize within for easy plotting
mutate(oxy_ct = oxy - mean(oxy)) %>%
ungroup()
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> pivot_longer: reorganized (Median, 1st, 3rd) into (series, oxy) [was 27x4, now 81x3]
#> ungroup: no grouping variables
#> group_by: one grouping variable (series)
#> mutate (grouped): new variable 'oxy_ct' (double) with 81 unique values and 0% NA
#> ungroup: no grouping variables
# Find which depths to calculate average oxygen on
wm_depth %>% filter(series == "1st") %>% summarise(mean_depth = mean(depth))
#> filter: removed 54 rows (67%), 27 rows remaining
#> summarise: now one row and one column, ungrouped
wm_depth %>% filter(series == "3rd") %>% summarise(mean_depth = mean(depth))
#> filter: removed 54 rows (67%), 27 rows remaining
#> summarise: now one row and one column, ungrouped
annual_oxy <- pred_grid %>%
drop_na(oxy) %>%
filter(depth > 29 & depth < 61) %>% ###
#filter(depth > 40 & depth < 60) %>% ###
group_by(year) %>%
summarize(median_oxy = median(oxy)) %>%
ungroup() %>%
rename("oxy" = "median_oxy") %>%
mutate(series = "Median (env.)",
oxy_ct = oxy - mean(oxy))
#> drop_na: no rows removed
#> filter: removed 187,542 rows (73%), 71,037 rows remaining
#> group_by: one grouping variable (year)
#> summarize: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> rename: renamed one variable (oxy)
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#> new variable 'oxy_ct' (double) with 27 unique values and 0% NA
oxygen_series <- bind_rows(wm_oxy, annual_oxy)
plot_w_oxy <-
ggplot(oxygen_series, aes(year, oxy, color = series, group = series, fill = series, linetype = series)) +
stat_smooth(method = "gam", formula = y ~ s(x, k = 4), se = FALSE, size = 0.8) +
geom_point(size = 1.3, alpha = 0.8) +
scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
scale_linetype_manual(values = c(1, 1, 1, 2)) +
scale_color_manual(values = pal, name = "Quantiles") +
scale_fill_manual(values = pal) +
guides(fill = "none", linetype = "none", color = guide_legend(nrow = 2, title.position = "top", title.hjust = 0)) +
labs(y = expression(paste("O" [2], " [ml/L]", sep = "")), x = "Year",
color = "") +
theme_plot(base_size = 11) +
theme(legend.position = c(0.25, 0.12),
legend.key.size = unit(0.25, "cm"),
legend.background = element_rect(fill = NA),
plot.margin = unit(c(0, 0, 0, 0), "cm"))
plot_w_oxy
# zoom in..
plot_w_oxy +
scale_y_continuous(breaks = seq(0, 10, by = 0.1))
#coord_cartesian(ylim = c(4.5, 5))
#
oxygen_series %>% as.data.frame()
## Temperature
# Plot temperature in 2006
temp_plot <- plot_map_raster +
geom_raster(data = filter(predict_mcod, year == "2006"), aes(x = X * 1000, y = Y * 1000, fill = temp)) +
scale_fill_viridis(option = "B", direction = -1) +
labs(fill = "Temperature [°C]", y = "") +
theme_plot(base_size = 11) +
theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.key.height = unit(0.1, "cm"),
legend.key.width = unit(0.35, "cm"),
legend.position = c(0.8, .085),
legend.direction = "horizontal",
legend.background = element_rect(fill = NA),
legend.text = element_text(size = 5),
legend.title = element_text(size = 6)) +
guides(fill = guide_colorbar(title.position = "top", title.hjust = 1)) +
geom_sf(size = 0.1) +
annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = lab_size) +
annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = lab_size) +
annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = lab_size) +
annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = lab_size) +
annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = lab_size, angle = 75) +
annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = lab_size, angle = 75)
#> filter: removed 249,002 rows (96%), 9,577 rows remaining
# Plot biomass-weighted temperature
# Weighted total
wm_temp <- predict_mcod %>%
group_by(year) %>%
summarise("Median" = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.5)),
"1st" = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.25)),
"3rd" = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.75))) %>%
pivot_longer(cols = c("Median", "1st", "3rd"),
names_to = "series", values_to = "temp") %>%
ungroup() %>%
group_by(series) %>% # standardize within for easy plotting
mutate(temp_ct = temp - mean(temp)) %>%
ungroup()
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> pivot_longer: reorganized (Median, 1st, 3rd) into (series, temp) [was 27x4, now 81x3]
#> ungroup: no grouping variables
#> group_by: one grouping variable (series)
#> mutate (grouped): new variable 'temp_ct' (double) with 81 unique values and 0% NA
#> ungroup: no grouping variables
# Find which depths to calculate average oxygen on
wm_depth %>% filter(series == "1st") %>% summarise(mean_depth = mean(depth))
#> filter: removed 54 rows (67%), 27 rows remaining
#> summarise: now one row and one column, ungrouped
wm_depth %>% filter(series == "3rd") %>% summarise(mean_depth = mean(depth))
#> filter: removed 54 rows (67%), 27 rows remaining
#> summarise: now one row and one column, ungrouped
annual_temp <- pred_grid %>%
drop_na(temp) %>%
filter(depth > 29 & depth < 61) %>% ###
#filter(depth > 40 & depth < 60) %>% ###
group_by(year) %>%
summarize(median_temp = median(temp)) %>%
ungroup() %>%
rename("temp" = "median_temp") %>%
mutate(series = "Median (env.)",
temp_ct = temp - mean(temp))
#> drop_na: no rows removed
#> filter: removed 187,542 rows (73%), 71,037 rows remaining
#> group_by: one grouping variable (year)
#> summarize: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> rename: renamed one variable (temp)
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#> new variable 'temp_ct' (double) with 27 unique values and 0% NA
temp_series <- bind_rows(wm_temp, annual_temp)
plot_w_temp <-
ggplot(temp_series, aes(year, temp, color = series, group = series, fill = series, linetype = series)) +
stat_smooth(method = "gam", formula = y ~ s(x, k = 4), se = FALSE, size = 0.8) +
geom_point(size = 1.3, alpha = 0.8) +
scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
scale_linetype_manual(values = c(1,1,1,2)) +
scale_color_manual(values = pal, name = "Quantiles") +
scale_fill_manual(values = pal) +
guides(fill = "none", linetype = "none", color = "none") +
labs(y = "Temperature [°C]", x = "Year", color = "") +
theme_plot(base_size = 11) +
theme(plot.margin = unit(c(0, 0, 0.9, 0), "cm"))
plot_w_temp
## Make combined plot
(bathy_plot | oxy_plot | temp_plot) / (plot_w_dep | plot_w_oxy | plot_w_temp) +
plot_annotation(tag_levels = 'A') &
theme(axis.text = element_text(size = 7),
axis.title = element_text(size = 7),
legend.text = element_text(size = 6),
legend.title = element_text(size = 6))
l <- (bathy_plot / oxy_plot / temp_plot) & theme(legend.text = element_text(size = 6),
legend.title = element_text(size = 7))
r <- (plot_w_dep / plot_w_oxy / plot_w_temp) & theme(legend.text = element_text(size = 8),
legend.title = element_text(size = 9))
(l | r) + plot_annotation(tag_levels = 'A')
#ggsave("figures/weighted_quantiles.pdf", width = 17, height = 11, units = "cm")
ggsave("figures/weighted_quantiles.pdf", width = 17, height = 20, units = "cm")
# Oxygen
wm_oxy_sd <- predict_mcod %>%
filter(depth > 29 & depth < 61) %>%
group_by(year, sub_div) %>%
summarise("weighted median" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.5)),
"un-weighted median" = median(oxy)) %>%
pivot_longer(c("weighted median", "un-weighted median"))
#> filter: removed 187,542 rows (73%), 71,037 rows remaining
#> group_by: 2 grouping variables (year, sub_div)
#> summarise: now 135 rows and 4 columns, one group variable remaining (year)
#> pivot_longer: reorganized (weighted median, un-weighted median) into (name, value) [was 135x4, now 270x4]
ggplot(wm_oxy_sd, aes(year, value, color = factor(name), fill = factor(name))) +
stat_smooth(method = "gam", formula = y ~ s(x, k = 4), size = 0.8) +
geom_point(size = 1.5) +
scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
scale_color_brewer(palette = "Accent", name = "") +
scale_fill_brewer(palette = "Accent") +
guides(fill = "none") +
labs(y = expression(paste("O" [2], " [ml/L]", sep = "")), x = "Year") +
theme_plot() +
facet_wrap(~sub_div) +
theme(aspect.ratio = 1,
legend.position = "bottom",
legend.key.size = unit(0.6, "cm"),
legend.text = element_text(size = 8),
plot.margin = unit(c(0, 0, 0, 0.2), "cm"))
ggsave("figures/supp/density/weighted_oxygen_subdiv.png", width = 6.5, height = 6.5, dpi = 600)
# Calculate weighted mean and quantiles for sd 25 only
predict_mcod %>%
filter(sub_div == 25 & year >= 2015) %>%
summarise("Median" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.5)),
"1st quartile" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.25)),
"3rd quartile" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.75))) %>%
as.data.frame()
#> filter: removed 245,134 rows (95%), 13,445 rows remaining
#> summarise: now one row and 3 columns, ungrouped
# Now compare delta change in oxygen with delta change in condition
# Depth
wm_dep_sd <- predict_mcod %>%
filter(depth > 29 & depth < 61) %>%
group_by(year, sub_div) %>%
summarise("weighted median" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.5)),
"un-weighted median" = median(depth)) %>%
pivot_longer(c("weighted median", "un-weighted median"))
#> filter: removed 187,542 rows (73%), 71,037 rows remaining
#> group_by: 2 grouping variables (year, sub_div)
#> summarise: now 135 rows and 4 columns, one group variable remaining (year)
#> pivot_longer: reorganized (weighted median, un-weighted median) into (name, value) [was 135x4, now 270x4]
ggplot(wm_dep_sd, aes(year, value, color = factor(name), fill = factor(name))) +
stat_smooth(data = filter(wm_dep_sd, name == "weighted median"), method = "gam", formula = y ~ s(x, k = 4), size = 0.8) +
geom_point(size = 1.5) +
scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
scale_color_brewer(palette = "Accent", name = "") +
scale_fill_brewer(palette = "Accent") +
guides(fill = "none") +
labs(y = "Depth [m]", x = "Year", color = "") +
theme_plot() +
facet_wrap(~sub_div) +
theme(aspect.ratio = 1,
legend.position = "bottom",
legend.key.size = unit(0.6, "cm"),
legend.text = element_text(size = 8),
plot.margin = unit(c(0, 0, 0, 0.2), "cm"))
#> filter (grouped): removed 135 rows (50%), 135 rows remaining
ggsave("figures/supp/density/weighted_depth_subdiv.png", width = 6.5, height = 6.5, dpi = 600)
# Now plot together
p1 <- ggplot(filter(wm_dep_sd, name == "weighted median"),
aes(year, value, color = factor(sub_div))) +
stat_smooth(data = filter(wm_dep_sd, name == "weighted median"), method = "gam", formula = y ~ s(x, k = 4), size = 0.8, se = FALSE) +
geom_point(size = 1.5) +
scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
scale_color_brewer(palette = "Dark2", name = "") +
guides(fill = "none") +
labs(y = "Depth [m]", x = "Year", color = "") +
theme_plot() +
theme(aspect.ratio = 1,
legend.position = "bottom",
legend.key.size = unit(0.6, "cm"),
legend.text = element_text(size = 8),
plot.margin = unit(c(0, 0, 0, 0.2), "cm"))
#> filter (grouped): removed 135 rows (50%), 135 rows remaining
#> filter (grouped): removed 135 rows (50%), 135 rows remaining
wm_oxy_sd2 <- filter(wm_oxy_sd, name == "weighted median")
#> filter (grouped): removed 135 rows (50%), 135 rows remaining
m1 <- lm(value ~ year*as.factor(sub_div), data = wm_oxy_sd2)
m2 <- lm(value ~ year+as.factor(sub_div), data = wm_oxy_sd2)
p2 <- ggplot(wm_oxy_sd2, aes(year, value, color = factor(sub_div))) +
stat_smooth(data = filter(wm_oxy_sd, name == "weighted median"), method = "gam", formula = y ~ s(x, k = 4), size = 0.8, se = FALSE) +
geom_point(size = 1.5) +
scale_x_continuous(breaks = c(1994, 2000, 2006, 2012, 2018)) +
scale_color_brewer(palette = "Dark2", name = "") +
guides(fill = "none") +
labs(y = "Depth [m]", x = "Year", color = "") +
theme_plot() +
theme(aspect.ratio = 1,
legend.position = "bottom",
legend.key.size = unit(0.6, "cm"),
legend.text = element_text(size = 8),
plot.margin = unit(c(0, 0, 0, 0.2), "cm"))
#> filter (grouped): removed 135 rows (50%), 135 rows remaining
p2
# Temporal trends of sprat by sd
pred_grid %>%
ggplot(aes(x = year, y = biomass_spr_sd, color = factor(sub_div))) +
stat_smooth(method = "gam", formula = y ~ s(x, k = 4), size = 1) +
geom_point(size = 2.5, alpha = 0.8) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2", name = "ICES subdivision") +
guides(fill = "none", linetype = "none") +
labs(y = "Sprat biomass [tonnes]", x = "Year",
color = "") +
theme_plot() +
theme(aspect.ratio = 1,
legend.position = "bottom",
legend.key.size = unit(0.6, "cm"),
legend.text = element_text(size = 8),
plot.margin = unit(c(0, 0.4, 0, 0), "cm"))
#> Warning: Removed 4693 rows containing non-finite values (`stat_smooth()`).
#> Warning: Removed 4693 rows containing missing values (`geom_point()`).
ggsave("figures/supp/sprat_subdiv.png", width = 6.5, height = 6.5, dpi = 600)
#> Warning: Removed 4693 rows containing non-finite values (`stat_smooth()`).
#> Removed 4693 rows containing missing values (`geom_point()`).
# Prepare prediction data frame
d2 <- d %>% drop_na(oxy)
nd_oxy <- data.frame(oxy = seq(min(d2$oxy), max(d2$oxy), length.out = 100))
nd_oxy <- nd_oxy %>%
mutate(year = 2003L,
depth_sc = 0,
oxy_sc = (oxy - mean(oxy))/sd(oxy),
temp_sc = 0)
# Predict
p_cond_oxy <- predict(m, newdata = nd_oxy, se_fit = TRUE, re_form = NA)
cond_oxy <- ggplot(p_cond_oxy, aes(oxy, exp(est),
ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se))) +
geom_line() +
geom_ribbon(alpha = 0.4) +
coord_cartesian(ylim = c(30, 220)) +
theme(aspect.ratio = 1) +
labs(x = expression(Oxygen[haul]))
# Prepare prediction data frame
nd_dep <- data.frame(depth = seq(min(d2$depth), max(d2$depth), length.out = 100))
nd_dep <- nd_dep %>%
mutate(year = 2003L,
depth_sc = (depth - mean(depth))/sd(depth),
oxy_sc = 0,
temp_sc = 0)
# Predict
p_cond_dep <- predict(m, newdata = nd_dep, se_fit = TRUE, re_form = NA)
cond_depth <- ggplot(p_cond_dep, aes(depth, exp(est),
ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se))) +
geom_line() +
geom_ribbon(alpha = 0.4) +
coord_cartesian(ylim = c(30, 220)) +
theme(aspect.ratio = 1) +
labs(x = expression(Depth[haul]))
# Prepare prediction data frame
nd_temp <- data.frame(temp = seq(min(d$temp), max(d$temp), length.out = 100))
nd_temp <- nd_temp %>%
mutate(year = 2003L,
depth_sc = 0,
oxy_sc = 0,
temp_sc = (temp - mean(temp))/sd(temp))
# Predict
p_cond_temp <- predict(m, newdata = nd_temp, se_fit = TRUE, re_form = NA)
cond_temp <- ggplot(p_cond_temp, aes(temp, exp(est),
ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se))) +
geom_line() +
geom_ribbon(alpha = 0.4) +
coord_cartesian(ylim = c(30, 220)) +
theme(aspect.ratio = 1) +
labs(x = expression(Temp[haul]))
p_cond_dep2 <- p_cond_dep %>%
mutate(variable = "Depth") %>%
dplyr::select(est, est_se, depth_sc, variable) %>%
rename(value = depth_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_cond_oxy2 <- p_cond_oxy %>%
mutate(variable = "Oxygen") %>%
dplyr::select(est, est_se, oxy_sc, variable) %>%
rename(value = oxy_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_cond_tem2 <- p_cond_temp %>%
mutate(variable = "Temperature") %>%
dplyr::select(est, est_se, temp_sc, variable) %>%
rename(value = temp_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
conds <- bind_rows(p_cond_dep2, p_cond_oxy2, p_cond_tem2)
ggplot(conds, aes(value, exp(est), ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
facet_wrap(~variable, ncol = 2) +
labs(y = "Biomass density",
x = "Scaled value") +
theme_plot()
ggsave("figures/supp/density/conditional_effects.png", width = 6.5, height = 6.5, dpi = 600)
knitr::knit_exit()